******** Set up ********
clear
cd "/Users/Jacob/Box Sync/Capital-punishment-religion-prejudice"

use "Data/gss-cappun-CLEAN-restrictVars.dta"

*ssc install asdoc
*ssc install outreg2

set scheme s2mono

graph set window fontface "Times New Roman"

******** Add labels ********
label variable punish2_dp "Support for death penalty"
label variable reltrad "Religious tradition"
label variable rr_index_new "Racial resentment"
label variable bible_literal "Biblical literalism"
label variable church_attend "Religiosity"
label variable partyid "Party identification"
label variable black "Black"
label variable other_race "Other race"
label variable coninc "Family income"
label variable south "South"
label variable age "Age"
label variable bach_degree "Bachelor's degree"
label variable female "Female"

*label define reltrad1 1 "Non-affiliate" 2 "Other religion" 3 "Jewish" 4 "Roman Catholic" 5 "Black Protestant" 6 "Mainline Protestant" 7 "Evangelical Protestant"
*label values reltrad reltrad1

label define rr_index_new1 0 "Low" 1 "Mod. low" 2 "Mod. high" 3 "High"
label values rr_index_new rr_index_new1

******** Clean variables ********
replace partyid = 3 if partyid == 7
gen log_coninc = ln(coninc)

******** Drop missing values ********
drop if punish2_dp == .
drop if reltrad == .
drop if rr_index_new == .
drop if bible_literal == .
drop if church_attend == .
drop if partyid == .
drop if black == .
drop if other_race == .
drop if log_coninc == .
drop if south == .
drop if age == .
drop if bach_degree == .
drop if female == .

******** Drop Black Protestants and Black respondents ********
drop if reltrad == 5
drop if black == 1

******** Factor analysis ********
alpha rr2_lackwill rr3_discrimination rr_single
factor rr2_lackwill rr3_discrimination rr_single

******** Descriptive statistics table ********
asdoc sum punish2_dp i.reltrad i.rr_index_new bible_literal church_attend partyid other_race log_coninc south age bach_degree female, save(Final-tables/summary-stats.doc) replace


******** RR dependent variable ********
* Model 1

asdoc poisson rr_index_new i.reltrad, save(Final-tables/rr-dv.doc) label nest replace cnames(Model 1)

* Model 2
asdoc poisson rr_index_new i.reltrad church_attend partyid other_race log_coninc south age bach_degree female i.year, save(Final-tables/rr-dv.doc) label nest cnames(Model 2)


******** Main results ********

* Model 1
asdoc logit punish2_dp i.reltrad bible_literal church_attend partyid other_race log_coninc south age bach_degree female i.year, or save(Final-tables/main-results.doc) label nest replace cnames(Model 1)


* Model 2
asdoc logit punish2_dp i.reltrad i.rr_index_new bible_literal church_attend partyid other_race log_coninc south age bach_degree female i.year, or label nest cnames(Model 2)


* Model 2
asdoc logit punish2_dp i.reltrad##i.rr_index_new bible_literal church_attend partyid other_race log_coninc south age bach_degree female i.year, or label nest cnames(Model 3)


******** UPDATED Predicted probability bar plots ********
logit punish2_dp i.reltrad##i.rr_index_new bible_literal church_attend partyid other_race log_coninc south age bach_degree female i.year, or
	  
estimate store a1 /** store in working memory **/
	  
	margins rr_index_new, post at(reltrad=(1))
	
	marginsplot, recast (bar) x(rr_index_new) plotopts(barw(.5)) ///
	ytitle("Predicted % support for CP") title("Non-affiliate", color(black)) name(nones, replace) ///
	xtitle("") ///
	xlabel(, nogrid) ///
	ylabel(.5 "50%" .6 "60%" .7 "70%" .8 "80%" .9 "90%", nogrid) ///
	plotregion(margin(b=1 l=3 t=0) style(none)) legend(off) ///
	graphregion(fcolor(white) lcolor(white)) ///
	plotregion(margin(b=1 l=3 t=3) style(none))
	
estimates restore a1

	  
	margins rr_index_new, post at(reltrad=(3))
	
	marginsplot, recast (bar) x(rr_index_new) plotopts(barw(.5)) ///
	ytitle("") title("Catholic", color(black)) name(cath, replace) ///
	xtitle("") ///
	ylabel(.5 "50%" .6 "60%" .7 "70%" .8 "80%" .9 "90%", nogrid) ///
	plotregion(margin(b=1 l=3 t=0) style(none)) legend(off) ///
	graphregion(fcolor(white) lcolor(white)) ///
	plotregion(margin(b=1 l=3 t=3) style(none))
	
estimates restore a1


	margins rr_index_new, post at(reltrad=(4))
	
	marginsplot, recast (bar) x(rr_index_new) plotopts(barw(.5)) ///
	ytitle("") title("Non-denominational Protestant", color(black)) name(non_denom, replace) ///
	xtitle("") ///
	xlabel(, nogrid) ///
	ylabel(.5 "50%" .6 "60%" .7 "70%" .8 "80%" .9 "90%", nogrid) ///
	plotregion(margin(b=1 l=3 t=0) style(none)) legend(off) ///
	graphregion(fcolor(white) lcolor(white)) ///
	plotregion(margin(b=1 l=3 t=3) style(none))
	
estimates restore a1


	margins rr_index_new, post at(reltrad=(6))
	
	marginsplot, recast (bar) x(rr_index_new) plotopts(barw(.5)) ///
	ytitle("") title("Mainline Protestant", color(black)) name(mp, replace) ///
	xtitle("") ///
	xlabel(, nogrid) ///
	ylabel(.5 "50%" .6 "60%" .7 "70%" .8 "80%" .9 "90%", nogrid) ///
	plotregion(margin(b=1 l=3 t=0) style(none)) legend(off) ///
	graphregion(fcolor(white) lcolor(white)) ///
	plotregion(margin(b=1 l=3 t=3) style(none))
	
estimates restore a1


	margins rr_index_new, post at(reltrad=(7))
	
	marginsplot, recast (bar) x(rr_index_new) plotopts(barw(.5)) ///
	ytitle("") title("Evangelical Protestant", color(black)) name(ep, replace) ///
	xtitle("") ///
	xlabel(, nogrid) ///
	ylabel(.5 "50%" .6 "60%" .7 "70%" .8 "80%" .9 "90%", nogrid) ///
	plotregion(margin(b=1 l=3 t=0) style(none)) legend(off) ///
	graphregion(fcolor(white) lcolor(white)) ///
	plotregion(margin(b=1 l=3 t=3) style(none))
	
	
graph combine nones cath non_denom mp ep, col(3) xsize(8) row(1) ycommon name(final_plot, replace) ///
	plotregion(margin(b=1 l=3 t=0) style(none)) ///
	graphregion(fcolor(white) lcolor(white)) ///
	plotregion(margin(b=1 l=3 t=3) style(none))

graph export "Final-figures/marginal-effects-final-noblack.pdf", replace


******** Coefficient equality tests ********
logit punish2_dp i.reltrad##i.rr_index_new bible_literal church_attend partyid other_race log_coninc south age bach_degree female i.year, or

margins rr_index_new, post at(reltrad=(1 2 3 4 6 7))

* Evangelical - Mainline (low RR) *
lincom _b[7bn._at#0.rr_index_new] - _b[6bn._at#0bn.rr_index_new]

* Evangelical - Catholic (low RR) *
lincom _b[7bn._at#0.rr_index_new] - _b[3bn._at#0bn.rr_index_new]

* Evangelical - non denom (low RR) *
lincom _b[7bn._at#0.rr_index_new] - _b[4bn._at#0bn.rr_index_new]

* Evangelical - BP (low RR) *
lincom _b[7bn._at#0.rr_index_new] - _b[5bn._at#0bn.rr_index_new]


* Evangelical - Mainline (high RR) *
lincom _b[7bn._at#3.rr_index_new] - _b[6bn._at#3bn.rr_index_new]

* Evangelical - Catholic (high RR) *
lincom _b[7bn._at#3.rr_index_new] - _b[3bn._at#3bn.rr_index_new]

* Evangelical - non denom (high RR) *
lincom _b[7bn._at#3.rr_index_new] - _b[4bn._at#3bn.rr_index_new]

* Evangelical - BP (high RR) *
lincom _b[7bn._at#3.rr_index_new] - _b[5bn._at#3bn.rr_index_new]



******** Mediation tests ********


**** Mediation tests (should be run between models 1 and 2)
recode reltrad 1=0 2=1 3/7=0, gen(other_trad)
recode reltrad 1/2=0 3=1 4/7=0, gen(r_catholic)
recode reltrad 1/3=0 4=1 5/7=0, gen(non_denom)
recode reltrad 1/4=0 5=1 6/7=0, gen(b_protestant)
recode reltrad 1/5=0 6=1 7=0, gen(main_protestant)
recode reltrad 1/6=0 7=1, gen(con_protestant)


// Other //

medeff (regress rr_index_new other_trad r_catholic non_denom b_protestant main_protestant con_protestant church_attend partyid black other_race log_coninc south age bach_degree female year) (logit punish2_dp other_trad r_catholic non_denom b_protestant main_protestant con_protestant rr_index_new  church_attend partyid black other_race log_coninc south age bach_degree female year), treat (other_trad) mediate (rr_index_new) sims (1000)

/* use % of total effect mediated*/
 *9.3% of the relationship is mediated*


// Catholic //

medeff (regress rr_index_new other_trad r_catholic non_denom b_protestant main_protestant con_protestant church_attend partyid black other_race log_coninc south age bach_degree female year) (logit punish2_dp other_trad r_catholic non_denom b_protestant main_protestant con_protestant rr_index_new  church_attend partyid black other_race log_coninc south age bach_degree female year), treat (r_catholic) mediate (rr_index_new) sims (1000)

/* use % of total effect mediated*/
 *34.7% of the relationship is mediated*
 
 
 // Non-denom Protestant //

medeff (regress rr_index_new other_trad r_catholic non_denom b_protestant main_protestant con_protestant church_attend partyid black other_race log_coninc south age bach_degree female year) (logit punish2_dp other_trad r_catholic non_denom b_protestant main_protestant con_protestant rr_index_new  church_attend partyid black other_race log_coninc south age bach_degree female year), treat (non_denom) mediate (rr_index_new) sims (1000)

/* use % of total effect mediated*/
 *17.5% of the relationship is mediated*

// Black Protestant //

medeff (regress rr_index_new other_trad r_catholic non_denom b_protestant main_protestant con_protestant church_attend partyid black other_race log_coninc south age bach_degree female year) (logit punish2_dp other_trad r_catholic non_denom b_protestant main_protestant con_protestant rr_index_new  church_attend partyid black other_race log_coninc south age bach_degree female year), treat (b_protestant) mediate (rr_index_new) sims (1000)

/* use % of total effect mediated*/
 *13.9% of the relationship is mediated*
 
 // Mainline Protestant //
 
 medeff (regress rr_index_new other_trad r_catholic non_denom b_protestant main_protestant con_protestant church_attend partyid black other_race log_coninc south age bach_degree female year) (logit punish2_dp other_trad r_catholic non_denom b_protestant main_protestant con_protestant rr_index_new  church_attend partyid black other_race log_coninc south age bach_degree female year), treat (main_protestant) mediate (rr_index_new) sims (1000)
 
/* Use % of total effect mediated*/
 *14.8% of the relationship is mediated*
 
 // Conservative Protestant //
  medeff (regress rr_index_new other_trad r_catholic non_denom b_protestant main_protestant con_protestant church_attend partyid black other_race log_coninc south age bach_degree female year) (logit punish2_dp other_trad r_catholic non_denom b_protestant main_protestant con_protestant rr_index_new  church_attend partyid black other_race log_coninc south age bach_degree female year), treat (con_protestant) mediate (rr_index_new) sims (1000)
 
/* Use % of total effect mediated*/
 *15.5% of the relationship is mediated*




